### Code to make figures for Fowler and Hall
### NOTE: this file should only be run *after* the Stata code
### the Stata code creates data files this code requires

## adjust working directory here
setwd("~/Dropbox/Fowler-Hall/Sharks/FowlerHall_Sharks_ReplicationData")

### install this package if necessary
### install.packages("readstata13")
library(readstata13)

### override default mgp in par to make better ticks
par.gh <- function(..., mgp=c(0,.05,.05)) {
  par(..., mgp=mgp)
}

### function to produce nice plots
### preserves most defaults through "..."
### defaults labels to blank, axes to missing, makes points look nice
plot.gh <- function(..., xlab="", ylab="", xaxt="n", yaxt="n", pch=21, bg="gray80", col="black",bor_col="white", x_seq = NULL,  y_seq = NULL, bty=NULL, rect_col="gray90", do_axes=NULL, do_seq=T, outer.box=F, outer.box.col="black", grid_col="white", point.size=1.3, remove.axis=F) {
  plot(..., xlab=xlab, ylab=ylab, xaxt=xaxt, yaxt=yaxt, pch=pch, bg=bg,  bty="n")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = rect_col, border= bor_col)
  # Draw sequences in background
  if(do_seq!=F){
    grid(col=grid_col, lty=2)
  }
  # Add axis
  if (!remove.axis) {
    axis(side=1,lwd=0, lwd.ticks=1, tck=-0.01, col="white", col.axis="gray20", cex.axis=.8)
    axis(side=2,lwd=0, lwd.ticks=1, tck=-0.01, las=1, col="white", col.axis="gray20", cex.axis=.8)
  }
  # Add points
  points(...,pch=21, bg=bg, col=col, lwd=1.5, cex=point.size)
  # Y/N add box around plot
  if(outer.box==T){
    box(col= outer.box.col)
  }
}

### Michael Gill's brilliant function to make transparent colors
makeTransparent<-function(someColor, alpha)
{
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}


#### make plots

data <- read.dta13("NJ_CountyData_RGraph.dta")
data.m <- data[data$machine==1,]

#### county level plot

pdf(file="county_graph.pdf", width=11, height=5)
par.gh(mfrow=c(1,2), oma=c(3,3,3,3), mar=c(1,1,1,3))
plot.gh(x=data$wilson1912, y=data$wilson1916, ylim=c(0.3,.7), xlim=c(0.3,0.7), col="gray90", bg="gray90")
cols <- ifelse(data$beach==1, "red", "gray50")
cols[data$coastal==1 & data$beach==0] <- "royalblue"
pts <- ifelse(data$beach==1, 15, 16)
pts[data$coastal==1 & data$beach==0] <- 17
points(x=data$wilson1912, y=data$wilson1916, pch=pts, col=cols, cex=1, lwd=1.2)
legend("topleft", bty="n", pch=c(15,17,16,1), col=c("red", "royalblue", "gray50", "black"), c("Beach Counties", "Other Coastal Counties", "Interior Counties", "Machine Counties"), cex=.8, pt.cex=c(1,1,1,2))
mtext(side=2, line=2, "Wilson 2-Party Vote Share, 1916", cex=.8)
mtext(side=3, line=.1, "Raw Data", font=2)
mtext(side=1, line=2, "Wilson 3-Party Vote Share, 1912", cex=.8)
points(x=data.m$wilson1912, y=data.m$wilson1916, cex=3)
plot.gh(x=data$wilson1912, y=data$wilson1916, ylim=c(-0.2,.2), xlim=c(-0.2,0.2), col="gray90", bg="gray90")
points(x=data$wilson1912res, y=data$wilson1916_nobeach, bg="gray70", col=cols, pch=pts, cex=1, lwd=1.2)
points(x=data$wilson1912res, y=data$wilson1916_beach, bg="red", col="red", cex=1, lwd=1.2, pch=15)
axis(side=1, lwd=0, cex.axis=.8)
mtext(side=3, line=.1, "Controlling for Machine Counties", font=2)
mtext(side=2, line=2, "Residualized Wilson Vote, 1916", cex=.8)
mtext(side=1, line=2, "Residualized Wilson Vote, 1912", cex=.8)
legend("topleft", bty="n", pch=c(15,17,16), col=c("red", "royalblue", "gray50"), c("Beach Counties", "Other Coastal Counties", "Interior Counties"), cex=.8)
dev.off()

### town-level plot

pdf(file="town_graph.pdf", width=9, height=5)
data <- read.dta13("NJ_TownData_RGraph.dta")
data$mat <- 0
data$mat[grep("Matawan", data$town)] <- 1
data.o <- data[data$county=="Ocean",]
cexs.o <- 1.2*(data.o$total1916 / max(data.o$total1916)) + .5
cols.o <- ifelse(data.o$beach==1, "red", "gray50")
cols.o[data.o$nearbeach==1] <- "royalblue"
pchs.o <- ifelse(data.o$beach==1, 15,16)
pchs.o[data.o$nearbeach==1] <- 17
par.gh(mfrow=c(1,2), oma=c(3,3,3,3), mar=c(1,1,1,1))
plot.gh(x=data.o$vote1912, y=data.o$vote1916,bg="gray90", xlim=c(.1,.8), ylim=c(.1,.7), col="gray90")
points(x=data.o$vote1912, y=data.o$vote1916, col=cols.o, pch=pchs.o, cex=cexs.o)
points(x=data.o$vote1912[data.o$town=="Sea Side Park"], y=data.o$vote1916[data.o$town=="Sea Side Park"], cex=2)
legend("topleft", bty="n", pch=c(15,17,16), col=c("red", "dodgerblue", "gray50"), c("Beach Towns", "Near-Beach Towns", "Non-Beach Towns"), cex=.8)
arrows(x0=.36, x1=.41, y0=.54, y1=.56, length=.1, angle=20)
mtext(side=2, line=2, "Wilson 2-Party Vote Share, 1916", cex=.9)
mtext(side=3, line=.1, "Ocean County, NJ", font=2)
text(x=.27, y=.54, "Sea Side Park", cex=.7)
cexs <- 1.2*(data$total1916 / max(data$total1916)) + .5
pchs <- ifelse(data$beach==1, 15,16)
pchs[data$nearbeach==1] <- 17
cols <- ifelse(data$beach==1, "red", "gray50")
cols[data$nearbeach==1] <- "royalblue"
plot.gh(x=data$vote1912, y=data$vote1916,bg="gray90", xlim=c(.1,.8), ylim=c(.1,.7), remove.axis=T, col="gray90")
points(x=data$vote1912, y=data$vote1916, col=cols, pch=pchs, cex=cexs)
points(x=data[data$mat==1,]$vote1912, y=data[data$mat==1,]$vote1916, cex=2)
axis(side=1, lwd=0, cex.axis=.8)
mtext(side=1, outer=T, line=1, "Wilson 3-Party Vote Share, 1912", cex=.9)
mtext(side=3, line=.1, "All NJ Counties", font=2)
legend("topleft", bty="n", pch=c(15,17,16,1), col=c("red", "dodgerblue", "gray50", "black"), c("Beach Towns", "Near-Beach Towns", "Non-Beach Towns", "Matawan Attacks"), pt.cex=c(1,1,1,2), cex=.8)
dev.off()





